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Abstract 

A methodology for solving two-point boundary value problems in phase space for Hamil- 
tonian systems is presented. Using Hamilton- Jacobi theory in conjunction with the canonical 
transformation induced by the phase flow, we show that the generating functions for this trans- 
formation solve any two-point boundary value problem in phase space. Properties of the gener- 
ating functions are exposed, we especially emphasize multiple solutions, singularities, relations 
with the state transition matrix and symmetries. Then, we show that using Hamilton's prin- 
cipal function we are also able to solve two-point boundary value problems, nevertheless both 
methodologies have fundamental differences that we explore. Finally, we present some applica- 
tions of this theory. Using the generating functions for the phase flow canonical transformation 
we are able to solve the optimal control problem (without an initial guess) , to study phase space 
structures in Hamiltonian dynamical systems (periodic orbits, equilibrium points) and classical 
targeting problems (this last topic finds its applications in the design of spacecraft formation 
trajectories, reconfiguration, formation keeping, etc.). 

1 Introduction 

One of the most famous two-point boundary value problems in astrodynamics is Lambert's 
problem, which consists of finding a trajectory in the two-body problem which goes through two 
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given points in a given lapse of time. Even though the two-body problem is integrable, no analytical 
solution has been found to this problem so far, and solving Lambert's problem still requires one to 
solve Kepler's equation, which has motivated many papers since 1650 5 . For a general Hamiltonian 
dynamical system, a two-point boundary value problem is solved using shooting methods combined 
with Newton iteration. Though very systematic, this technique requires a "good" initial guess for 
convergence and is not appropriate when several boundary value problems need to be solved. In 
order to design a change of configuration of a formation of n spacecraft, n\ two-point boundary value 
problems need to be solved ^H] , hence for a large collection of spacecraft the shooting method is not 
efficient. In this paper we address a technique which allows us to solve m boundary value problems 
at the cost of m function evaluations once generating functions for the canonical transformation 
induced by the phase flow are known. These generating functions are solutions of the Hamilton- 
Jacobi equation and for a certain class of problem they can be found offline, that is during mission 
planning. Moreover, the theory we expose allows us to formally solve any kind of two-point boundary 
value problem, that is, given a n-dimensional Hamiltonian system and 2n coordinates among the 
An defining two points in the phase space, we find the other 2n coordinates. The Lambert problem 
is a particular case of this problem where the dynamics is Keplerian, the position of two points are 
given and the corresponding momenta need to be found. Another instance of such a problem is the 
search for trajectories which go through two given points in the momentum space (i.e., the conjugate 
of the Lambert problem). Properties of the solutions found are studied, such as multiple solutions, 
symmetries and relation to the state transition matrix for linear systems. Then, we expose another 
method to solve two-point boundary value problems based on Hamilton's principal function and 
study how it compares to generating functions. Finally, we present direct applications of this theory 
through the optimal control problem and the study of some Hamiltonian dynamical systems. Solving 
the optimal control problem using generating functions was first introduced by Scheeres et al. |17j . 
we will review their method in this paper and expand it to more general optimal control problems. 
Applications to Hamiltonian dynamical systems were first studied by Guibout and Scheeres [§llll)| 
for spacecraft formation flight design and for the computation of periodic orbits. 



2 Solving a two-point boundary value problem 

In this section, we recall the principle of least action for Hamiltonian systems and derive the 
Hamilton- Jacobi equation. Local existence of generating functions is proved. We underline that we 
do not study global properties. In general, we do not know a priori if the generating functions will 
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be defined for all time and in most of the cases we found that they develop singularities. We refer 
the reader to [Tll^ir7ll51Pll3lll4| for more details on local Hamilton-Jacobi theory, [11121113] for global 
theory and (HIEE! and section l2*.3. 4l of this paper for a study of singularities. 

2.1 The Hamilton-Jacobi theory 

Let (P,uj,Xff) be a Hamiltonian system with n degrees of freedom, and H : P x R — > M the 
Hamiltonian function. In the extended phase space Pxl, we consider an integral curve of the 
vector field Xh connecting the points (qo,po,to) and (qi,pi,ti). The principle of least action reads: 

Theorem 2.1. (The principle of least action in phase space) The integral L pdq — Hdt has 

an extremal in the class of curve 7 whose ends lie in the n- dimensional subspaces (t = to,q = qa) 
and (t = ti, q = qi) of extended phase space. 

Proof. We proceed to the computation of the variation. 

f f ( dH dH \ 
$ / (PQ ~ H)dt = / ( qSp + pSq — Sq — 5p ) dt 

J-y J-y\ Oq ' 



- \pSq]l+ I 
J ~i 



dp 

W)T) — n -I 



dt (2.1) 



Therefore, since the variation vanishes at the end points, the integral curves of the Hamiltonian 
vector field are the only extremals. □ 

Remark 2.1. The condition for a curve 7 to be an extremal of a functional does not depend on the 
choice of coordinate system, therefore the principle of least action is coordinate invariant. 

Now let (Pi,wi) and (P2,0J2) be symplectic manifolds, 

Definition 2.1. A smooth map / : Pi x R — > P2 x R is a canonical transformation if and only if 

(1) - / is a C°°-diffeomorphism, 

(2) - / preserves the time, i.e., there exists a function g t such that f(x,t) = (gt(x),t), 

(3) - for each t, gt : Pi — > P2 as defined above is a symplectic diffeomorphism and / preserves the 
canonical form of Hamilton's equations. 

All three points in this definition are not independent but wc mention them for sake of clarity. 
It can be proved that if gt is symplectic then / is a diffeomorphism. Moreover, the third point 
of the definition differs from book to book. We chose Abraham's definition pQ but very often the 
third item reduces to "/ preserves Hamilton's equations" (Goldstein [7], Greenwood Arnold [2] 
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argues that this definition differs from the original definition, the third item should actually be u gt 
is symplectic" which implies, but is not equivalent to, "/ preserves the canonical form of Hamilton's 
equations" . 

Consider a canonical transformation / : (qi,Pi,t) i— > (Qi,Pi,t). Since Hamilton's equations are 
preserved, we have: 

W ' " dP l 

(2.2) 

• = _dK_ 

where K = K(Q. P, t) is the Hamiltonian of the system in the new set of coordinates. 

On the other hand, we have seen that the principle of least action is coordinate invariant. Hence: 

S f 1 (X>&-#(«,p,t) j dt = (2.3) 



to 



\i=l 



5 J \T^PiQi-K{Q,P,t)jdt = Q (2.4) 

From Eqns. 12.31 - l2~4l we conclude that the integrands of the two integrals differ at most by a 
total time derivative of an arbitrary function F: 

n n 

Pidqt -Hdt = Y^ P J d Qj ~ Kdt + dF ( 2 - 5 ) 

i=l 3 =1 

Such a function is called a generating function for the canonical transformation / and is, a priori, 
a function of both the old and the new variables and time. The two sets of coordinates being con- 
nected by the 2n equations, namely, f(qi, ■ ■■ , q n ,Pi, ■ ■ • ,p n ,t) = (Qi, •■• , Qn, Pi,- ■ , Pn,t), F can 
be reduced to a function of 2n+l variables among the 4n+l. Hence, we can define 4™ generating func- 
tions that have n variables in Pi and n in P2. Among these are the four kinds defined by Goldstein \7\, 
Fi(qi,--- ,q n ,Qi,"' ,Qn,t), F 2 (qi,--- ,q n ,Pi,"' ,Pn,t), F 3 (p 1} --- ,p n , Qi, ■ • ■ ,Q n ,t) and 
Fa(pi, ■ ■ ■ ,Pn,Pi, ■ ■ ■ ,P n ,t). 

Let us first consider the generating function Fi(q, Q, t). The total time derivative of F\ reads: 

dF!(q,Q,t) = ^Idqi + f^^-dQi + ^-dt (2.6) 

■ - , oq % ^— ' dQi at 

i=i ^ 3=1 ^ 
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Hence Eq. 12. 51 yields: 

E [Pi ~ ~ Hdt = E (p >' + m-) d ®i ~ Kdt + d ~w dt ( 2J ) 

i—l 1 j— l ^ 

Assume that (q, Q 1 1) is a set of independent variables, then Eq. 12.71 is equivalent to: 

P* = -^(g.Q,*) (2-8) 

oqi 

If (g, Q) is not a set of independent variables, we say that F\ is singular. 

Now let us consider more general generating functions. Let (ij., • ■ • , i p )(ip+i, • • • , i n ) and 
(&!,•■• , fc r )(fc r+ i, • • • , k n ) be two partitions of the set (1, •• • ,n) into two non-intersecting parts such 
that ii < • • • < i p , ip+1 < • • • < i n , k± < ■ ■ ■ < k r and k r+ i < ■ ■ ■ < k n and define I p = («i, • • • , i p ), 

4 = (ip+l)" - >*n)) #r = (fel,-" and AT,. = (fer+l,'-' If 

(<ll P ,PT p >QK r ,PK r ) = fen - ' ' >5i„)Pi P+ i)- ' ' ,Pi n ,Qki,' ' ' , Qfc r , -Pfc r+ i , ' ' ' ,-PfeJ 

are independent variables, then we can define the generating function Fj p K r '- 
F i P ,K r (<li P , Pl p ,QK r ,PR r , t) = F iln,- ■ ■ ,Qi p ,Pi p +i,-- >Pin>Qk!,--- ,Qk r ,Pk r+1 r ■ ■ ,Pk n ,t) (2.11) 

Expanding dFj p ,K r yields: 

dFj K r , dFj k , dFj k,. ,„ dFj x r JT , dFj K 

a— 1 a a— p+1 a— 1 a a— r+1 a 

(2.12) 

and rewriting Eq. 12.51 as a function of the linearly independent variables leads to: 

p n r n 

^Pijli a - £ «*- d P*« -Hdt = Y^ PkJQk a - £ QkJPk a - Kdt + dF Ip , Kr (2.13) 

a— 1 a— p+1 a— 1 a— r+1 

where Fj p ,K r — Fi + Sa=r+i Qk a Pk a ~ Y^l= P +i liaPia This last relation defines the Legendre trans- 
formation, which allows one to transform one generating function into another. 
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Then Eq. |2~T31 reads: 



a=l ^ ka a=r+l ka 



dt 



a=l 



a=p+l 



(2.14) 



which is equivalent to: 

dF 

pi, = -^f^(qi P ,Pi p ,QK r ,P Rr ,t) (2.15) 

OF 

% = ^-(qi p ,p Ip ,QK r ,PK r ,t) (2.16) 

OF 

p K r = --~%^{qi f ,Pi,,QK r ,PR r ,t) (2.17) 

OF 

Qk t = -g^(qi P ,Pi p ,QK r ,P Rr ,t) (2.18) 

T ^,^ dFi K r dFj k dFi k t dFi K r dFj k 

For the case where the partitions are (1, • • • , n)() and ()(1, • • • , n) (i.e.,p = n and r — 0), we 
recover the generating function F2, which verifies the following equations: 

F)F 

Pl = ^(q,P,t) (2.20) 
oqi 

dF 

Qi = gp;(<l,P,t) (2-21) 

dF 9 dF? OF? 

K{ -l,P, t) . H( 9 ,-^, t) + - s l (2.22, 



The case p — and r — n corresponds to a generating function of the third kind, F3: 



<1< 
Pi 



FlF 



r) F 

OPi 

F)F 

8F 3 A ,dF 3 



(2.23) 
(2.24) 
(2.25) 
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Finally, if p = and r — 0, we obtain F^: 

% = ^(p,P,t) (2.26) 

op 

•'•<> - "<f- (, + f < 2 - 28 > 

2.2 The phase flow is a canonical transformation 

In the following we focus on a specific canonical transformation, the one induced by the phase 
flow. Let <I> t be the flow of an Hamiltonian system: 

$ t :P P 

(<lo,Po) >-> (®l{qo,Po)=q(qo,Po,t),®%(qo,Po)=p{qo,Po,t)) (2.29) 
Then, the phase flow induces a transformation <fi on P x R defined as follows: 

0:(g O ,2>o,i)^($t(?o,Po),i) (2.30) 
Theorem 2.2. TTie transformation cj) induced by the phase flow is canonical. 

Proof. The proof of this theorem can be found in Arnold [2], it is based on the integral invariant of 

Poincare-Cartan. □ 

For such a transformation, (Q,P) represents the initial conditions of the system (qo,po), the 
Hamiltonian function if is a constant that can be chosen to be and the equations verified by the 
generating function Fi v .k t become: 

pip = g ^ r {qi P , pr p , qo Kr ,po Rr ,t), t) (2.3i) 

qi p = — Q^-iqip'Pip'qoK^Poz^t),?) (2.32) 

Po Kr = — g q P o - (qip,PT p ,qo Kr ,POg r ,t),t) (2.33) 

OFj k 

Qo Rr = ^ '- - (qi p , Pi p ,Qa Kr , Po Rr ,t), t) (2.34) 

tt, dFl v ,K r 9Fr t K r ,s OFj Kr ,„ „„ N 
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The last equation is often referred to as the Hamilton- Jacobi equation. To solve this equation, one 
needs boundary conditions. At the initial time, position and momentum (g,p) are equal to the initial 
conditions ((foiPo)- Hence, Fj pt K r must generate the identity transformation at the initial time. 

2.3 Properties of the canonical transformation induced by the phase flow 

In this section we study the properties of generating functions for the phase flow canonical 
transformation. First we show that they solve a two-point boundary value problem, and then we 
prove a few results on singularities, symmetries and differentiability. In particular, we relate the 
generating functions and the state transition matrix for a linear system. 

2.3.1 Solving a two-point boundary value problem 

Consider two points in phase space, Xq = (qa,pa) and X\ = (q,p), and two partitions of (1, • • • , n) 
into two non-intersecting parts, (ii,-- - , i p ) (ip+i , • • • ,i n ) and (fci,-- - , fc r )(fc r +i, • • ■ , k n ). A two- 
point boundary value problem is formulated as follows: 

Given 2n coordinates (q^,-" , ?i fl R t+ i, ■ • ■ ,Pi n ) and (q 0ki . ■ ■ ■ , qo kr ,Po kr+1 , ■ - ■ ,Po kn ), find the re- 
maining 2n variables such that a particle starting at Xq will reach X\ in T units of time. 

From the relationship defined by Eqns. l2~T5l ETlTfl |2~T7I and ETTHl we see that the generating 
function Fi pi k t solves this problem. Lambert's problem is a particular case of boundary value 
problem where the partitions of (1, • • • ,n) are (1, • • • ,n)() and (1, • • • ,«)()■ Though, given two 
positions <?/ and go and a transfer time T, the corresponding momentum vectors are found from the 
relationships verified by F% : 

Pi = -7^—(q,qo,T) 

oqi 

Po t = -^(q^o^) (2.36) 

2.3.2 Existence and properties of the generating functions 

In the first section we proved the existence of a generating function using the assumption that 
its variables are linearly independent. This is not always true at every instant. As an example let 
us look at the harmonic oscillator. The equations of motion are given by: 



q(t) = qocos(ujt) + po I 'u> sin(u>t) 
p(t) = — qoLO sin(ujt) +pocos(ujt) 



(2.37) 
(2.38) 



9 



At T = 2ir/ui + 2kir, we have q(T) = q a , that is (q,qo) are not independent variables and the 
generating function Fi is undefined at this instant. We say that Fi is singular at T. We now prove 
that at least one of the generating functions is not singular at every instant. 

Proposition 2.3. Consider the flow $t of an Hamiltonian system tfi : (qo,Po,t) i— ► (^t(Qo>Po)>t)i 
where 

$t : (3o, Po) ($t((7o,Po) = q(qo,Po,t),$t(<lo,Po) =p(qo,Po,t)) 
For every t, there exists two subsets of cardinal n of the set (1, • ■ • , 2ri), I n and K n , such that 

det(||) ^0 (2.39) 

where $ t (q,p,q ,p ) = {$l(q,qo,Po) = Q ~ ®l(qo,Po), $?(p, qo,Po) = P - (<7o,Po)) and z = (q ,Po) 
Proof. To prove this property, we only need to notice that & t is a diffeomorphism, i.e., $t i€J is 
an injection, therefore, there exists at least one n-dimensional subspace on which the restriction of 
<&t ieIn is a diffeomorphism. □ 

Theorem 2.4. At every instant, at least one generating function is well-defined. Moreover, when 
they exist, generating functions define local C°° -diffeomorphism. 

Proof. From the previous theorem, there exists I n and K n such that dct ( ) ^ 0. 

Without any loss of generality and for simplicity, let the partition be /„ = (1, • • • , n) = K n , then we 
have: 

^{ ^dqT^ ) * ° (2 - 40) 



Moreover, <t] verifies: 



$ 1 t (q,q , Po ) = Q (2.41) 



From the local inversion theorem there exists a local diffeomorphism fi in a neighborhood of (qo,Po) 
such that qo = fi(q,Po)- In addition, the flow defines p as a function of (qo,Po), i- e -j replacing q by 
fi(q,Po), we obtain (q ,p) = (fi(q,Po), h{q,Pa)) where f 2 (q,Po) = ®1{fi{q,Po),Po)- This equation 
is equivalent to the two equations verified by F 2 , hence /1 = and f 2 = This proves that 

F 2 exists and since $ t is C°° , F 2 defines a local C°°-diffeomorphism from (q,po) to (p, qo). □ 

Remark 2.2. The theorem above can be stated for generating functions associated with an arbitrary 
canonical transformation, not only the one induced by the phase flow. To proceed the above proof 
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we only required that the flow defines a C^-diffeomorphism, this property is shared by all canonical 
transformations. 

Through the harmonic oscillator example, we saw that a generating function may become singu- 
lar. We now characterize singularities and give a physical interpretation to them. 

Proposition 2.5. The generating function Fi v> k t is singular at time t if and only if 



From the above theorem, we deduce that a generating function is singular when there exists 
multiple solutions to the boundary value problem. In the harmonic oscillator example, whatever the 
initial momentum is, the initial position and position at time T = 2n / uj + 2k-K are equal. 

Finally, if the Hamiltonian function is independent of time, the system is reversible and therefore 
the generating functions Fi p .K r and Fx r ,i p are similar in the sense that there exists a diffcomorphism 
which transforms one into the other. In particular, they develop singularities at the same instant. 
If p = n and r = 0, we obtain that F 2 and F 3 are similar. 

2.3.3 Linear systems theory 

In this section we particularize the theory developed above to linear systems. The following 
developments have implication in the study of relative motion and in optimal control theory as 
we will see later. Further, using linear systems theory, we are able to characterize singularities of 
generating functions using the state transition matrix. 

Hamilton-Jacobi equation When studying the relative motion of two particles, one often lin- 
earizes the dynamics about the trajectory of one of the particles (called the nominal trajectory) and 
then uses a linear approximation of the dynamics to study the motion of the other particle relative 
to the nominal trajectory (perturbed trajectory). Thus, the study of relative motion reduces to the 
study of a time-dependent linear Hamiltonian system, i.e., a system with a quadratic Hamiltonian 




(2.42) 



where I = {i G I p } [J{n + i,i G Ip) and J = {j G K r ) [J{n + j, j G K r ). 



Proof. The proof proceeds as the previous one, it is also based on local inversion theorem. 



□ 
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function without any linear terms UJ: 



1 yT I H qq {t) H qp (t) 

>Hpq(t) H pp (t) ^ 



(2.43) 



where Xh = f Ap ) ^ s ^ ne reia ti ve state vector. Guibout and Scheeres ^ proved that the generating 
functions for the phase flow transformation must then be quadratic without any linear terms, that 
is, if we take F2 for example: 



F 2 = -Y 1 



K Fl(t) F^{t) / 



Y 



(2.44) 



where Y = ^ £ p q Q ^ and ^ \ is the relative state vector at initial time. We also point out that 
both matrices defining Hh and F2 are symmetric. Then Eq. 12.201 reads: 



Ap 



dF 2 
dAq 



(2.45) 



Substituting into Eq. 12. 2 21 yields 



,1. 



r / 



Y T < 



FUt) T 




I 

FUt) T 



H qq (t) H qp (t) 

Hpqit) H pp (t) 



I 

Fh(t) FUt) 



Y (2.46) 



Though the above equations have been derived using F 2 , they are also valid for Fx (replacing 
Y = (a/ ) b y y = ( Ago )) smce ^1 anc ^ F 2 solve the same Hamilton- Jacobi equation (Eqns. |2.10| 

1 For the canonical transformation induced by the phase flow, we have seen that K = 
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I2.22|) . Eauation l2.4t)l is equivalent to the following 4 matrix equations: 

Fff(t) + H qq {t) + H qp {t)Fli 2 {t) + Fli 2 {t)H pq {t) + Fli 2 {t)H pp {t)Fli 2 {t) = 



Fli{t) + H qp {t)Fl 2 2 {t) + Fli 2 (t)H pp (t)F^ 2 (t) = 

(2.47) 

Fk\t) + F 2 \ 2 {t)H pq {t) + F^ 2 {t)H pp {t)Fli 2 {t) = 



Fl 2 2 {t) + F 2 \ 2 {t)H pp {t)Fl 2 2 {t) = 

where we replaced F^ by F^ to signify that these equations are valid for both F\ and F2 and recall 
that F 2 i = F\ 2 . A similar set of equations can be derived for any generating function Fi pt K r , 
here we only give the equations verified by F% and -F4: 

F?f(t) + H pp (t) - H pq (t)Fl { \t) - Fl { \t)H qp {t) + Fli\t)H qq {t)Flf{t) = 



Fh\t) - H pq {t)Fli\t) + Fti\t)H qq (i)Fl 2 \t) = 

(2.48) 

^ 2 3 i 4 « ~ F^(t)H qp (t) + Fli\t)H qq (t)F?f(t) = 

F^{t) + Fk\t) (t)H qq (t)F? 2 4 (t) = 

The first equations of Eans 12.471 and 171351 are Ricatti equations, the second and third are non- 
homogeneous, time varying, linear equations and are equivalent to each other (i.e., transform into 
each other under transpose), and the last are just a quadrature. 

Perturbation matrices Another approach to the study of relative motion of spacecraft is to use 
the state transition matrix. This method was developed by Battin 0] for the case of a spacecraft 
moving in a point mass gravity field. Let <& be the state transition matrix which describes the 
relative motion: 

' i9 U( AW | (2.49, 
Ap / I A Po 



where $ 



$ qp 



<5>pq $pp 
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From the state transition matrix, Battin 0] defines the fundamental perturbation matrices C 
and C as: 



(7 = $ $ _1 

VP qp 



(2.50) 



That is, given Ap — 0, CAq = Ap and given A<7 = 0, CAq — Ap. He shows that for relative 
motion of a spacecraft in a point mass gravity field these matrices verify a Ricatti equation and are 
therefore symmetric. Using the generating functions for the canonical transformation induced by 
the phase flow, wc immediately recover these properties and also show that they are verified for any 
relative motion of two particles in a Hamiltonian dynamical system. 
From Eqns. |2~27)I and |2~2"T1 



Ap = 



dF 2 
dA Po 

F^Aq 

dF 2 
dAq 
F^Aq- 



F? 2 A Po 



Fi 2 A Po 



(2.51) 



(2.52) 



Solving for (Aq, Ap) yields: 



Aq = Fl x Aq - F% x Fi 2 Ap 

Ap = F^Fi^Aqo + (ff 2 - F^F^ 1 F 2 2 2 )Ap 



(2.53) 
(2.54) 



From the above equations we are able to link the state transition matrix to the generating function 
F 2 . 



<s> qp 


= -Flii 


%q 


- ^21 


%p 


= F? 2 -F, 


% q 


= Fl.Fl- 



We conclude that 



C = $ p9 $-i = p» 



(2.55) 
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In the same manner, but using Pi, we can show that: 

C = = F\ x (2.56) 

Thus, we have shown that C and C are symmetric by nature (as P^ 2 is symmetric by definition) , 
and moreover that they verify the Ricatti equation given in Eq. 12.471 

Singularities of generating functions and state transition matrix From Eqns. I2.2UI 12.211 

and 12.441 we derive a relationship between terms of F 2 and some coefficients of the state transition 
matrix: 



Thus: 



Ap 



dF 2 
dAp 



Ap 



Aq 



Aq a = 



F^Aq- 


■ F? 2 A Po 




but we 


also have 




%q^qq 


Ag + ($ pp - <f> qp )Ap 




dF 2 
dAp 
F^Aq- 


f P 2 2 2 Ap 


(2.57) 


but we 


also have 








(2.58) 


*?i = 


Vqq 1 


(2.59) 


Fh = 


(f> — (J) <f)~ let) 

^PP ^pq^qq ^qp 


(2.60) 


F 2 \ = 


$-1 

qq 


(2.61) 


F 22 




(2.62) 



(2.63) 



We conclude that if Q qq is singular, F 2 is also singular. The same analysis can be achieved for 
the other generating functions, and we find that: 

• Pi is singular when <J>„„ is singular, 
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• i<3 is singular when <fr pp is singular, 

• F4 is singular when <J> pg is singular. 

These results can be extended to other generating functions, but requires us to work with another 
block decomposition of the state transition matrix. 

2.3.4 On singularities of generating functions 

We have proved local existence of generating functions and mentioned that they may not be glob- 
ally defined for all time. Using linear system theory we were able to predict where the singularities 
are and to interpret their meaning as multiple solutions to the two-point boundary value problem. 
In this section we extend our study to singularities of nonlinear systems. 

Lagrangian submanifold Consider an arbitrary generating function Fj k t - Then Eqns. 12. 151 
12.181 define a 2n-dimensional submanifold called a canonical relation U£ of the 4n-dimensional 
symplectic space P\ x P 2 . I n addition, since the new variables (Q,P) (or (qo,Po)) do n °t appear 
in the Hamilton-Jacobi equation 12.351 we may consider them as parameters. In that case Eqns. 
12.151 and 12.161 define an n-dimensional submanifold of the symplectic space Pi called a Lagrangian 
submanifold |19|. The study of singularities can be achieved using either canonical relations 1 or 
Lagrangian submanifolds piUl], 

Theorem 2.6. The generating function Fi vi k t is singular if and only if the local projection of the 
canonical relation C defined by Eqns. \2.15H2.T8i onto {gi p ,Pj v1 Qk t , P{i r ) is n °t a l° ca l diffeomor- 
phism. 

Moreover, the projection of such a singular point is called a caustic. If one works with Lagrangian 
submanifolds then the previous theorem becomes 

Theorem 2.7. The generating function 2 Fj p ^K r is singular if the local projection of the Lagrangian 
submanifold defined by Eqns. \2.15\ and \2.1b\ onto (qi p ,pr p ) is not a local diffeomorphism. 

In light of these previous theorems, we can give a geometrical interpretation to theorem 12.41 on 
the existence of generating functions. Given a canonical relation C (or a Lagrangian submanifold) 
defined by a canonical transformation, there exists a 2n-dimensional (or n-dimensional) submanifold 
M. of P\ x P 2 (or Pi) such that the local projection of C onto Ai is a local diffeomorphism. 

2 We consider here that the generating function is function of n variables only, and has n parameters. 



16 



Study of caustics To study caustics two approaches, at least, are possible depending on the 
problem. A good understanding of the physics may provide information very easily. For instance, 
consider the two body problem in dimension 2, and the problem of going from one point A to a point 
B, symmetric with respect to the central body, in a certain lapse of time. The trajectory that links 
A to B is an ellipse whose perigee and apogee are A and B. Therefore, there are two solutions to 
this problem depending upon which way the particle is going. In terms of generating functions, we 
deduce that F3 is nonsingular (there is a unique solution once the final momentum is given) but F% is 
singular (existence of two solutions) and the caustic is a fold 3 . The other method to study caustics 
consists in using a known nonsingular generating function to define the Lagrangian submanifold 
C and then study its projection. A very illustrative example is given by Ehlers and Newman 
they treat the evolution of an ensemble of free particles whose initial momentum distribution is 
P = \+qi using the Hamilton-Jacobi equation and generating functions for the phase flow canonical 
transformation. They are able to solve the problem analytically, that is, identify a time at which F± 
is singular, find the equations defining the Lagrangian submanifold using F3 and study its projection 
to eventually find two folds. Nevertheless, such an analysis is not always possible as solutions to 
the Hamilton-Jacobi equation are usually found numerically, not analytically. In the remainder of 
this section, we focus on a class of problem that can be solved numerically for which we are able to 
characterize the caustics. 

Suppose we are interested in the relative motion of a particle, called the deputy, whose coordinates 
are (q,p) with respect to another one, called the chief, whose coordinates are (q ,p ), both moving 
in an Hamiltonian field. If both particles stay "close" to each other, we can expand (q,p) as a Taylor 
series about the trajectory of the chief. The dynamics of the relative motion is described by the 
Hamiltonian function Hh [H|: 



p=2 1 2n ' ' ' ' OQnOPl • ' ' Opn 



Vlln=P 

(2.64) 

where X h = (Aq,Ap), Aq = q — q° and Ap = p — p°. Since Hh has infinitely many terms, we 
are usually not able to solve the Hamilton-Jacobi equation but we can approximate the dynamics 
by truncating the series Hh in order to only keep finitely many terms. Suppose N terms are kept, 
then we say that we describe the relative motion using an approximation of order N. Clearly, the 

3 Maps from R 2 into R 2 have two types of stable singularities, folds and cusps. However, only folds have two 
antecedents, cusps have three. 
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greater N is, the better our approximation is to the nonlinear motion of a particle about the nominal 
trajectory. When an approximation of order N is used, we look for a generating function Fi Pt K r as 
a polynomial of order N in its spatial variables with time dependent coefficients. The Hamilton- 
Jacobi equation reduces to a set of ordinary differential equations that we integrate numerically. 
Once F[ p ^K r is known, we find the other generating functions from the Legrendre transformation, 
at the cost of a series inversion. If a generating function is singular, the inversion does not have a 
unique solution, the number of solutions characterizes the caustic. To illustrate this method, let us 
consider the following example. 

Motion about the Libration point L 2 in the Hill three-body problem Consider a spacecraft 
moving about and staying close to the Libration point Li in the Hill three-body problem (See the 
appendix for a description of the Hill three-body problem). Its relative motion with respect to Li is 
described by the Hamiltonian function (Eq. I2.64[) and approximated at order N by truncation 
of terms of order greater than N in the Taylor series defining Hh . Using the algorithm developed by 
Guibout and Scheeres |S] we find the generating functions for the canonical transformation induced 
by the approximation of the phase flow, that is, the Taylor series expansion up to order N of the 
exact generating function about the Libration point Li. 

Fi{q x ,q y ,Po^Po y ,t) = fu(t)ql + fl 2 {t)q x q y + fi 3 (t)q x p 0x + ftS)q x p Qy 

fUtWy + fis(t)q y PoAt) + fUt)lyP% 

flS)pl x + fhitfPOvPay + fLi^POy + r (.Qx,q y ,Po^Pa v ,t) (2.65) 

where (q,p, qo,Po) are relative position and momenta of the spacecraft with respect to Li at t and 
at £0 ; the initial time, and r is a polynomial of degree N in its spatial variables with time dependent 
coefficients and without any quadratic terms. At T = 1.6822, F\ is singular but Fi is not. Eqns. 
|2~2TJI and |2~2T1 reads: 

p x = 2ft x {T)q x + fli{T)q y + fUT)p^ + fu(T)po y + F>ir(q x ,q y ,po x ,po y ,T) (2.66) 
Py = lli( T )lx + 2 f22( T )ly + fli( T )Po* + fli( T )P% + F>ir(q Xl q y ,p 0j; ,po y , T) (2.67) 
= f?3( T )^ + f23( T )<ly + 2 f33( T )PQ,+fu( T )Po v +F>3r(q x ,q y ,p 0j; ,po y ,T) (2.68) 
1o y - fu(T)q x + fii(T)q y + .f 3 2 4 (7> ,+2 f! i (T)p 0y + D 4 r(q x ,q y ,p 0x ,p 0y ,T) (2.69) 
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where Z^r represents the derivative of r with respect to its i th variable. Eqns. I2.66I2.6T?! define a 
canonical relation C. By assumption i*\ is singular, therefore the projection of C onto (g, go) is n °t 
a local diffeomorphism and there exists a caustic. The theory developed above provides a technique 
to study this caustic using F%. Eqns. I2.6612.6TH provide p and go as a function of (q,po), but to 
characterize the caustic we need p and po as a function of (g, go)- Fx being singular, there are 
multiple solutions to this problem, and one valuable piece of information is the number k of such 
solutions. To find p and po as a function of (g, go) we can first invert equations 12.681 and 12.691 to 
express po as a function of (g, go) and then plug this relation into Eqns. I2.66l and l2.67l The first step 
requires a series inversion that can be proceeded using the technique developed in |15| by Moulton. 
Let us rewrite Eqns. 12.681 and !2. 691 

2 /33( T )Po x + fL(T)p 0y = Qo x ~ fx 3 ( T )Qx - f$ 3 (T)q y - D 3 r(q x ,q y ,p 0a:7 p 0y ,T) (2.70) 
fli{ T )pQ* + 2/ 4 2 4 (7>o H = qo y - fi 4 (T)q x - /| 4 ( T )?s/ - D A r(q x ,q y ,p 0lc ,p Qy ,T) (2.71) 

The determinant of the coefficients of the linear terms on the left hand side is zero (otherwise there 
is a unique solution to the series inversion) but each of the coefficients is non zero, that is, we can 
solve for po x as a function of (po y , go x , qo y ) using eauation l2.70l Then we plug this solution into Eq. 
12.711 and we obtain an equation of the form 

R(Po y ,qo x ,qo y )=0 (2.72) 

that contains no terms in po y alone of the first degree. In addition, R contains a non zero term of 
the form otp\ y , where a is a real number. In this case, Weierstrass proved that there exist 2 solutions 
pl y and p\ y to Eq. 12.721 that is, the caustic is a fold. 

In the same way, we can study the singularity of F± at initial time. At t = 0, F2 generates the 
identity transformation, hence /33(C)) = /34(D) = /34(D) = /44(D) = 0. This time there is no nonzero 
first minor, and we find that there exists infinitely many solutions to the series inversion. Another 
way to see this is to use the Legendre transformation: 

F%(q,qo,t) = F 2 (q,p ,t) - q Q p Q (2.73) 

As t tends toward 0, (q,p) goes to (go,Po) and F2 converges toward the identity transformation 
qPo -^t^o goPo- Therefore, as t goes to 0, -Fi also goes to 0, i.e., the projection of C onto (g, go) 
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reduces to a point. 

There are many other results on caustics and Lagrangian submanifolds that go beyond the 
scope of this paper. Study of the Lagrangian submanifold at singularities is "the beginning of 
deep connections between symplectic geometry, geometric optics, partial differential equations, and 
Fourier integral operators." (R. Abraham Q]), we refer to Abraham [I] and references given therein 
for more information on this subject. Let us now come back to two-point boundary value problems. 

So far we have studied the generating functions associated with the canonical transformation 
induced by the phase flow and showed they formally solve any two-point boundary value problem. 
Nonetheless, for Hamiltonian dynamical systems there exists another function, called Hamilton's 
principal function, that solves the same problem and thus for completeness we discuss it. In this 
section we introduce this function and show how it compares to the generating functions for the 
canonical transformation induced by the phase flow. 

2.4 Hamilton's principal function 

Though generating functions are used in this paper to solve boundary value problems, they 
have been introduced by Jacobi and mostly used thereafter as fundamental functions which can 
yield all the equations of motion by simple differentiations and eliminations, without integration. 
Nevertheless, it was Hamilton who first hit upon the idea of finding such a fundamental function, 
he proved its existence in geometrical optics (i.e., for time independent Hamiltonian systems) in 
1834 and called it characteristic function JT]. The year later, he published a second essay [T2] 
on systems of attracting and repelling points in which he showed that the evolution of dynamical 
systems is characterized by a single function called Hamilton's principal function: "The former 
Essay contained a general method for reducing all the most important problems of dynamics to 
the study of one characteristic function, one central or radical relation. It was remarked at the 
close of that Essay, that many eliminations required by this method in its first conception, might 
be avoided by a general transformation, introducing the time explicitly into a part S of the whole 
characteristic function V ; and it is now proposed to fix the attention chiefly on this part S, and to 
call it the Principal Function." (William R. Hamilton, in the introductory remarks of "Second essay 
on a General Method in Dynamics" \1'2\) 



20 



2.4.1 Hamilton's principal function to describe the phase flow 

As with generating functions, Hamilton's principal function may be derived using the calculus 
of variations. Consider the extended action integral: 

A= [ 1 ipq' + Pt t')dT (2.74) 



TO 



under the auxiliary condition K(q 1 t,p,pt) — 0, where q' = dq/dr, pt is the momentum associated 
with the generalized coordinates t and K = p t + H. 

Define a line element da for the extended configuration space (q, t) by 

da = Ldt = Lt'dr (2.75) 

Then, we can connect two points (go, t ) and (q\, t{) of the extended configuration space by a shortest 
line 7 and measure its length from: 

A = I da= I Lt'dr (2.76) 

J *y J ~y 

The distance we obtain is function of the coordinates of the end-points and is called Hamilton's 
principal function: W(qo,to, qi,t\). 

From calculus of variations we know that the variation of the action A can be expressed as 
a function of the boundary terms if we vary the limits of the integral: 

8 A = p Q 8q Q + p to 5t - p\8qi - p tl 5h (2.77) 

On the other hand we have: 

dW dW dW dW 
5 A = SW(q Q ,t , quh) = —6q Q + -^-5t + —5 gi + —6h (2.78) 

oq ot oqi dti 



that is: 



dW 

Po = -^—(qo,to,qi,ti) (2.79) 
oqo 

dW 

Pi = — 5— (go,*o, Qi, h) (2.80) 
oqi 



'Note that the geometry established by this line element is not Riemannian I13| 
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and 

dW dW 

— ( qo ,t ,qx,ti) + H{q ,—,t Q ) = (2.81) 
at oq 

dW dW 
— 3j-(9b. *o, 91, h) +H(q 1 ,- — ,t 1 ) = (2.82) 

where K has been replaced by pt + H. As with generating functions of the first kind, Hamilton's 
principal function solves boundary value problems of Lambert's type through Eqns. I2.79l and l2.80l 
To find W, however, we need to solve a system of two partial differential equations (Eqns. 12.811 and 

EH3- 

2.4.2 Hamilton's principal function and generating functions 

In this section we highlight the main differences between generating functions for the canonical 
transformation induced by the phase flow and Hamilton's principal function. For sake of simplicity 
we compare Fi(q,q ,t) and W(q, t, go, to). 

Calculus of variation Even if both functions are derived from calculus of variations, there are 
fundamental differences between them. To derive generating functions we used the principle of 
least action with the time t as independent variables whereas we increase the dimensionality of the 
system by adding the time t to the generalized coordinates to derive Hamilton's principal function. 
As a consequence, generating functions generates a transformation between two points in the phase 
space, i.e., they act without passage of time whereas Hamilton's principal function generates a 
transformation between two points in the extended phase space, i.e., between two points in the 
phase space with different times. This difference may be viewed as follows: Generating functions 
allow to characterize the phase flow given an initial time, to (i.e., to characterize all trajectories 
whose initial conditions are specified at to), whereas Hamilton's principal function does not impose 
any constraint on the initial time. The counterpart being that Hamilton's principal function must 
satisfy two partial differential equations (Eq. I2.81l defines W as a function of to and Eq. 12. 821 defines 
W as a function of t%) whereas generating functions satisfy only one. 

Moreover, to derive the generating functions fixed endpoints are imposed, that is we impose 
the trajectory in both sets of variables to verify the principle of least action. On the other hand, 
the variation used to derive Hamilton's principal function involves moving endpoints and an energy 
constraint. This difference may be interpreted as follows: Hamilton's principal function generates a 
transformation which maps a point of a given energy surface to another point on the same energy 
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surface and is not defined for points that do not lie on this surface. As a consequence of the energy 
constraint we have: 

d 2 W 

7jT7H = ° (2.83) 

As noticed by Lanczos ^3] > "this is a characteristic property of the W- function which has no equiv- 
alent in Jacobi's theory" . On the other hand, generating functions map any point of the phase space 
into another one, the only constraint is imposed through the principle of least action (or equivalcntly 
by the definition of canonical transformation): we impose the trajectory in both sets of coordinates 
to be Hamiltonian with Hamiltonian function H and K respectively. 

Fixed initial time In the derivation of Hamilton's principal function dto may be chosen to be 
zero, that is, the initial time is imposed. Hamilton's principal function loses its dependence with 
respect to to, Eq. 12.811 is trivially verified and Eq. 12.831 does not hold anymore, W and F\ become 
equivalent. 

Finally, in |12j Hamilton also derives another principal function Q(po,to,pi,ti) which compares 
to W as i*4 compares to Fx, the derivation being the same we will not go through it. 

To conclude, Hamilton's principal function appears to be more general than the generating 
functions for the canonical transformation induced by the phase flow. On the other hand, to solve 
a two-point boundary value problem, initial and final times are specified and therefore, any of these 
functions will identically solve the problem. To find Hamilton's principal function, we need to solve 
two partial differential equations whereas only one need to be solved to find the generating functions. 
For this reason, generating functions will be used in the following examples. 



3 Applications 

3.1 Solving the optimal control problem using the generating functions 

The use of the generating functions to solve an optimal control problem has first been addressed 
by Scheeres, Guibout and Park [17j . They suggested an indirect approach for evaluating the initial 
adjoints without initial guess. In the present paper, we review their approach and generalize it to a 
wider class of problem. 
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Problem formulation Assume a dynamical system described by: 

x = f(x,u,t) (3.1) 

x(t = 0) = x (3.2) 

where u is the control variable, and u E K m . An optimal control problem is formulated as 
follows: 

mmK(x(t f ))+ L(x,u,t)dt (3.3) 

where tf is the known final time. This formulation is called the Bolza formulation. Other formula- 
tions are possible and completely equivalent 

mmK(x(tf)) Mayer formulation (3.4) 

u 

C l S _ 

min / L(x,u,t)dt Lagrange formulation (3.5) 

u Jto 

Further, some final conditions may be specified. For instance, suppose that k final conditions 
are given for the final state, i.e., 

^j(x(t f ),t f ) = ;6(l-t) (3.6) 

Necessary conditions Define the Hamiltonian function H: 

H(x,p,u,t) = p T x + L(x,u,t) (3.7) 

where p G R™ is the costate vector. Applying the Pontryagin principle allows one to find the optimal 
control: 

u = argmini? (x,p, u, t) (3.8) 

u 

Then the necessary conditions for optimality are given by: 

x = ——(x,p,u,t) (3.9) 



dp 
81 

dx 



OH 

P = —z-(x,p,u,t) (3.10) 
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To integrate these 2n differential equations we need In boundary conditions: n + k are specified in 
the problem statement, the other n — k are given by the transversality conditions: 

P(*/)-^(*/) = " T ^(*(*/)) (3- 11 ) 
where v is a fc-dimensional vector. 

Solving the optimal control using the generating functions In the following, we are making 
two assumptions which may be relaxed in future research. 

1. One can solve for u as a function of (x,p) using Eq. 13.81 that is, we can define a new Hamiltonian 
function H(x,p,t) = H(x,p,u(x,p,t),t). 

2. One can eliminate the j/'s in Eq. 13.111 so that Eq. 13. Ill becomes 

Pi{tf)=Pf< VtG(fc,n) (3.12) 

and transform Eq. 13. 61 into: 

x j (tf)) = x fj j- i !•••/,;• (3.13) 

Then, solving the optimal control problem is equivalent to find the solutions (x,p) satisfying: 

OH 

i = -g^iX'P't) ( 3 - 14 ) 
OH 

p = - — (x,p,t) (3.15) 

with boundary conditions 

x{t = 0) = Xq 

Xi (tf) = x h Vie(l,---,fc) (3.16) 
Pi{tf)) = Pft Vie (fc,-- - ,n) 

These equations define a two-point boundary value problem and hence arc usually difficult to 
solve because they generally require an estimate of the initial (or final) state, which usually has no 
physical interpretation. An indirect approach can be developed to solve this problem, namely, the 
use of the generating function Fj n ,K k (xo 1 , ■ ■ ■ , xo n , Xf x , • ■ • ,%f h ,Pf h+1 , • • ■ >Pf n )- Eqns. 12.151 12T1UI 
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and 12.1 71 solves the boundary value problem and hence the optimal control problem: 



POi 

Pfi 



9F In , Kk 
dx 0i 

d Ph 
dF In , Kk 
dx fi 



(3.17) 
(3.18) 
(3.19) 



In the case where k = n, that is initial and final states of the system are specified, the generating 
function that must be used to solve the boundary value problem is F\. In that case, Park and 
Scheeres ^5] showed that F\ satisfies the Hamilton- Jacobi-Bellmann equation. 

Particular case: The linear quadratic problem Assume the dynamics of the system is linear: 



x(t) = A(t)x(t) + B(t)u(t) 



(3.20) 



and the cost function J is quadratic: 



J = ^[Mx(tf ) - m f ] T Q f [Mx(tf) - m f ] + 













u 1 


)t 








(3.21) 



and Q is symmetric positive semi-definite, R and Qf are symmetric positive definite. Moreover, 

( ' 
. . . x 

define L to be L = i ( X T u T 




Using previous notations, we define the Hamiltonian function H: 

H(x,p, u) — p T x + L(x. u) 

From eauation l3.8l we get 

u = -R- 1 B T p - R^^x 
Substituting u in Eqns. 13.91 and 13 . 1 Ul yields : 



(3.22) 



(3.23) 



H(x,p) = H(x,p, -R- 1 B T p - R^^x) 



(3.24) 
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and 



Ax + B(-R- 1 B T p - R^^x) 



(3.25) 



p = -{A T p + Qx + N{-R- T B T p-R- 1 N T x)) 



(3.26) 



Boundary conditions for this problem are still given by equations 13.161 Since the Hamiltonian 
function defining this system is quadratic, this problem is often solved using the state transition 
matrix. We have seen previously that, in linear systems theory, both generating functions and the 
state transition matrix are equivalent. Moreover, to compute the generating function or the state 
transition matrix, four matrix equations of dimension n need to be solved. Therefore, both methods 
are exactly equivalent for the linear quadratic problem. Finally, another method to solve the linear 
quadratic problem is to apply Ricatti transformation to reduce the problem to two matrix ordinary 
differential equations, a Ricatti equation and a time-varying linear equation. An analogy can be 
drawn between these two equations and the ones verified by the generating function. 

3.2 Finding periodic orbits using the generating functions 

Another application of the generating functions for the canonical transformation induced by the 
phase flow is to search for periodic orbits. This application was first presented by Guibout and 
Scheeres jlOj . we review their methodology in this paper and refer to |10| for more details and 
additional examples. 

3.2.1 The search for periodic orbits: a two-point boundary value problem 

The main idea is to transform the search for periodic orbits into a two-point boundary value 
problem that can be handled using generating functions. For a periodic orbit of period T, both 
position and momentum take the same values at t and at t + kT, k G Z. In terms of initial 
conditions, this reads: 



For a dynamical system with n degrees of freedom Eqns. 13.271 and 13.281 can be viewed as In 
equations of 2n + 1 variables, the initial conditions (qo,Po) and the period T. To solve such a 
problem, for each trial (qo,po,T) one needs to integrate the equations of motion and check if the 



q(T) = qo 



(3.27) 



P(T) = Po 



(3.28) 
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2n equations are verified, and if they are not try again. On the other hand, Eqns. 13.271 and 13.281 
can also be viewed as a two-point boundary value problem. Suppose the initial momentum p and 
the position at time T, q, are given, then Eqns. 13.271 and 13.281 define 2n equations with 2n + 1 
variables, the initial position go, the momentum at time T, p, and the time period T. Solutions to 
these equations characterize all periodic orbits. The idea now is to use the generating functions for 
the phase flow transformation to solve this problem. Depending on the two-point boundary value 
problem we choose to characterize periodic orbits, different generating functions can be used. In the 
following we will only deal with generating functions of the first and second kind, but this theory 
can be readily generalized to any kind of generating functions. 

3.2.2 Solving the two-point boundary value problem 

Generating functions of the first kind The generating function F\ allows us to solve a two- 
point boundary value problem for which initial position and position at time T are given. The 
solution to this problem is found using Eqns. 12.81 and [231 

P = ^j±(q t qo,T) (3.29) 

Po = -^-(q,qo,T) (3.30) 
oqo 

On the other hand, the boundary value problem that characterizes periodic orbits is defined by 
equations 13 .271 and 13 .281 Hence, combining these four equations yields: 

dF\ 

po = -^—(q = qo,qo,T) (3.31) 

oqo 

P(T) = -g 1 (q = qo,qo,T) (3.32) 

That is, since p{T) = po: 

dF 1 dFi 

-£—(9 = qo,qa,T) + -^—(q = qo,qo,T) = (3.33) 
dq Oqo 

Eq. 13.331 defines n equations with n + 1 variables, (qo,T), it is an under-determined system, and 
hence we often focus on one of the two following problems: 

1. Search in time domain: Given a point in the position space qo, find all periodic orbits going 
through this point and their associated momentum. Eq. 13.331 defines n equations of a single 
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variable T. Taking the norm of the left hand side yields: 

dF x dFi 

W-Q^il = 9o, qo,T) + -g-^{q = Qo,qo,T)\\ = (3.34) 

Eq. 13.341 is a single equation of one variable that can be solved graphically. To find the 
corresponding momentum, we can use either Eq. I2.2UI or Eq. 12.211 

dFi 

Po = — (q = qo,qo,T) (3.35) 
Oqo 

p = -5— {q = qo,qo,T) (3.36) 



Both equations provide the same momentum since Eq. 13. 341 is equivalent to \\p — po\\ =0 and 
is satisfied. 

2. Search in position space: Find all periodic orbits of a given period. Eq. 13.331 reduces to a 
system of n equations with n unknowns, qo. For dynamical systems with n degrees of freedom 
the solution may be represented on a n-dimensional plot. In practice, solving this problem 
graphically is efficient only for systems with at most 3 degrees of freedom. For Hamiltonian 
systems with more than 3 degrees of freedom, Newton iteration or an equivalent method can 
be used. When a solution to Eq. 13.331 is obtained, then we use Eq. 12.81 or 12.91 to find the 
corresponding momentum: 

dFi 

po = -^—(q = qo,qo,T) (3.37) 

Oqo 
dF 

p = -Q^-(q = qo,qo,T) (3.38) 

Generating function of the second kind The search for periodic orbits can also be solved 
using a generating function of the second kind. The main difference with the use of F% is that the 
system of equations we need to solve does not reduce to a system of n equations and n functions 
evaluations (we must solve 2n equations). 

The generating function F% allows us to solve a two-point boundary value problem for which the 
initial momentum and the position at time T are given. The solution to this problem is found using 
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Eqns. |2~2U1 and |2~2T1 

dF-2 

P = -^(g.po.T) (3.39) 
Qo = ^(g.Po.T) (3.40) 

On the other hand, the boundary value problem is defined by eauations l3.27l and l3.28l Combining 
these four equations yields: 



Po = p{T) 
dF 2 



(q,Po,T) (3.41) 



= ^l(q,p ,T) (3.42) 
opo 

The system of equations 13.411 and !3. 421 contains 2n equations with 2n + 1 variables, and therefore 
is under-determined. As with F\, we consider two main problems, we either set the time period or 
n coordinates of the point in the phase space. 

3.2.3 Examples 

To illustrate the theory developed above, let us consider the Hill three-body problem and let 
us find periodic orbits about the Libration point Li using the generating function of the first kind 
F\. To compute Fx, we use the algorithm developed by Guibout and Scheeres |Hj that computes 
the Taylor series expansion of the generating functions about a given trajectory, called the reference 
trajectory. In this example the reference trajectory is the equilibrium point L 2 and we compute the 
Taylor series up to order 6. Since we are working with series expansion, we will only find periodic 
orbits that stay within the radius of convergence of the series, not all periodic orbits. 

Search in time domain: Find all periodic orbits going through the point 5 (0.01,0). We have 
seen that this problem can be handled using Eq. 13. 341 which is one equation with one variable, T. In 
Figure^we have plotted the left-hand side of Eq. 13. 341 as a function of time, we obtain a continuous 
curve whose points have a particular significance. Let a; be a point on that curve whose coordinates 
are x — (t x ,Ap). The trajectory whose initial conditions are qo = (0-01,0), po — — j^-(qo,qo,t x ) 
comes back to its initial position after a time t x but the norm of the difference between its initial 



5 We use normalized units, for the Sun-Earth-spacecraft system 0.01 units of length represents about 21, 500fcm 
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momentum and its momentum at time t x is Ap. Hence, any point on the curve whose coordinates 
are (t x ,0) represents a periodic orbit (not only the trajectory comes back to its initial position at 
t x but the norm of the difference between the momenta at initial time and at t x is zero, i.e., the 
trajectory comes back to its initial state at t x ). In figure ^ we observe that there exists a periodic 
orbit of period T = 3.03353 going through the point (0.01, 0). The corresponding momenta is found 
using either Eq. TFM ov Eq. E^Tl and is p = p = (0, -0.0573157). 

Search in position space: Find all periodic orbits of period T = 3.0345. To solve this problem 
we use Eq. 13.331 which is a system of two equations with two variables (go x , qo y ) . In Fig. [5] wc 
plot solutions to each of these two equations and then superimpose them to find their intersection, 
which is the solution to Eq. 13.331 The solution is a closed curve, i.e., a periodic orbit of the given 
period. By plotting the solutions to Eq. 13.331 for different periods, we generate a map of a family 
of periodic orbits around the Libration point. In Figure |3| we plot the solutions to Eq. 13.331 for 
t = 3.033 + 0.0005n, «e{0, •••,9}. 

3.3 Study of equilibrium points 

The generating functions can also be used to study properties of equilibrium points of an Hamil- 
tonian dynamical system. First, we have proved the equivalence between the state transition matrix 
and the generating functions describing relative motion in linear system theory, therefore, linear 
terms in the Taylor series expansion of the generating functions about the equilibrium point provide 
information on the characteristic time and stability as does the state transition matrix. The other 
terms can be used to study the geometry of center, stable and unstable manifolds far from the 
equilibrium points where the linear approximation does not hold anymore (but within the radius of 
convergence of the Taylor series) . The study of center manifolds is a direct application of the previ- 
ous section as is readily seen from the example we provided. To find stable and unstable manifolds 
we propose a technique that uses generating functions to solve initial value problems, not two-point 
boundary value problems. Historically, generating functions were introduced by Jacobi and used 
thereafter to solve initial value problems, hence the following technique is not new. We mention it 
to show that one is able to fully describe an equilibrium point with only knowledge of the generating 
functions. 

The idea is to propagate the trajectory of a point that is "close" to the equilibrium point and 
on the linear approximation of the stable (unstable) manifold. Even though this method to find 
unstable and stable manifolds is not exact, it is fairly accurate and often used. We then reduce the 
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search for hyperbolic manifolds to an initial value problem that can be solved using any generating 
functions. For simplicity let us consider F 2 . At the linear level, a point on the unstable (stable) 
manifold has coordinates ((fthPo) = (ctu,a\u) where a <C 1, A is the characteristic exponent and u 
is the eigenvector defining the unstable (stable) manifold. Eq. 12.211 defines q(t) implicitly: 

dF 2 

au = q = -5 — (q,a\u,t) 

Once q(t) is found, we find p(t) from Eq. 12.201 As t varies, (q(t),p(t)) describes the hyperbolic 
manifolds. 

3.4 Design of spacecraft formation trajectories 

The last application we present concerns the design of a formation of spacecraft. This is again a 
direct application of the theory developed in this paper, first introduced by Guibout and Scheeres 
[S]. This application relies on the fact that the relative dynamics of two particles evolving in a 
Hamiltonian dynamical system is Hamiltonian, hence the Hamilton- Jacobi theory is applicable. To 
illustrate the use of generating functions, let us study an example. We consider a constellation of 
spacecraft located at the Libration point L 2 of Hill's three-body problem. At a later time t = tf, 
we want the spacecraft to lie on a circle surrounding the libration point at a distance of 108, OOOfcm. 
What initial velocity is required to transfer to this circle in time t / , and what will the final velocity 
be once we arrive at the circle? The answer will depend, of course, on where we arrive on the circle. 
In general, this problem must be solved repeatedly for each point on the circle we wish to transfer 
to and each transfer time. In our example we only need to compute the generating functions F± to 
be able to compute the answer as an analytic function of the final location. The method to solve 
this problem proceeds as follows: We first compute F\ then we compute the solution to the problem 
of transferring from L 2 to a point on the final circle where 2 parameters may vary, the transfer time 
and the location on the circle. Then we look at solutions which minimize the total fuel cost of the 
maneuver, that is, which minimize the sum of the norm of the initial momentum and the norm of 
the final momentum, ^/|Apo| 2 + |Ap| 2 . We assume zero momentum in the Hill's rotating frame at 
the beginning and end of the maneuver. While not a realistic maneuver, we can use it to exhibit 
the applicability of our approach. 

Figures [3] and show the value of ^/ 1 Apo 1 2 + |Ap| 2 as a function of position in the final 
formation at different times 6 . We notice three tendencies: 

^Define the final position of the spacecraft as Aq = Aqq where Aq = 108, OOOfcm and q is the unit vector pointing 
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1. For t less than the characteristic time, no matter which direction the spacecraft leaves L2, it 
costs essentially the same amount of fuel to reach the final position and stop (figure 0}. 

2. For t larger than the characteristic time, but less than 47 days the curve describing 
^/|Apo| 2 + |Ap| 2 shrinks along a direction 80° from the x-direction. Thus, placing a spacecraft 
on the final circle at an angle of 80° or 260° from the x-axis provides the lowest cost in fuel 
(figure 0) . 

3. For t larger than 47 days, the curve describing y^Apop + |App shrinks along a direction 
perpendicular to the previous one, at an angle of ~ 170° with the x-axis and expands along 
the 80° direction. Thus, there exists an epoch for which placing a spacecraft on the final circle 
at an angle of 170° or 350° from the x-axis provides the lowest final cost, this happens for 
t — 88 days (figures and EJ. 

To conclude, we see the optimal transfer time to the final circle changes as a function of location 
on the circle. While this is to be expected, our results provide direct solutions for this non-linear 
boundary value problem. 

We now make a few additional remarks to emphasize the advantage of our method. First, 
additional spacecraft do not require any additional computations. Hence, our method to design 
optimal reconfiguration is valid for infinitely many spacecraft in formation. Second, now that we 
have computed the generating functions around the libration point, we are able to analyze any 
reconfiguration around the libration point at the cost of evaluating a polynomial function 7 . Finally, 
if the formation of spacecraft is evolving around a base which is on a given trajectory, we can linearize 
about this trajectory, and then proceed as in the above examples to study the reconfiguration 
problem. 

Conclusions 

This paper describes a novel application of Hamilton-Jacobi theory. We are able to formally 
solve any nonlinear two-point boundary value problem using generating functions for the canonical 
transformation induced by the phase flow. Many applications of this method are possible, and we 
have introduced a few of them, and implemented them successfully. Nevertheless, the method we 

towards the location of the final circle. Then, figures 14161 represent iy|Apo| 2 + |Ap| 2 <j 

7 This is especially valuable for missions involving spacecraft that stay close to Li since the generating functions 
in this region can be computed during mission planning. Then any targeting problem or reconfiguration design can 
be achieved at the cost of a function evaluation 
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propose is based on our ability to obtain generating functions, that is to solve the Hamilton- Jacobi 
equation. In general such a solution cannot be found, but for a certain class of problem an algorithm 
has been developed 9 that converges locally in phase space. A typical use of this algorithm would 
be to study the optimal control problem about a known trajectory, to find families of periodic orbits 
about an equilibrium point or in the vicinity of another periodic orbit, and to study spacecraft 
formation trajectories. 



Appendix I: The circular restricted three-body problem and 
Hill's three-body problem 

The circular restricted three-body problem is a three-body problem where the second body is 
in circular orbit about the first one and the third body has negligible mass 0- The coordinate 
system is centered at the center of mass of the two bodies with mass and the Hamiltonian function 
describing the dynamics of the third body is: 



H(q x ,q v ,p x ,p v ) = \{p 2 x +pI) +p x q v - q x p v , 1 = = , M (3.43) 

1 ^/{qx + ^ 2 + q 2 y y \q x - i + m) 2 + ql 

where q x — x, q y — y, p x = x — y and p y — y + x. Equations of motion for the third body can be 
found from Hamilton's equations: 

= x - {1 - M) «* + w + ^ - " «& - 1 + w + w (3 - 44) 

" + 2X = y - {1 -"\(q X+ ,r + ql)^- f \(q X ~l + lr + q^ ^ 

There are five equilibrium points for this system, called the Libration points. L2 is the one whose 
coordinates are (1.01007,0) for a value of fx = 3.03591 ■ 10~ 6 . 

If the first body has a larger mass than the second one we can expand the equations of motion 
about /1 = 0. Then, shifting the coordinate system so that its center is the second body yields Hill's 
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formulation of the three-body problem. The equations are motion are: 



x-2y = + 3a; (3.46) 

y + 2x = -± (3.47) 



where r 2 = x 2 + y 2 . 

The Lagrangian then reads: 



Hence, 

dL 

Px = jp- 
oq x 



(3.48) 



1 13 

L{q,q,t) = - {ql+q 2 y ) + . + -q x - {q x q y - q y q x ) (3.49) 

kl + ql 1 



= 4x - q y (3.50) 
dL 
dq y 

= q y + q x (3.51) 



From Eqns. 13 . 491 f5TSUl and l3~5*Tl we obtain the Hamiltonian function H: 



H(q,p) = p x q x +p y q y -L 



M +P 2 y ) + (%Px ~ ixPy) ~ r= + ^ q v ~ 2q2 ^ ( 3 - 52 ) 

^ /"2 _i_ n 2 £ 



l,o , 1 I 

Iql + q' y 

There are two equilibrium points for this system, called libration points. Their coordinates are 
»V3 n \ „j r //nV> 



L!(-(i) 1/J ,0) and L 2 ((±) 1/a ,0) 



. lr 

0.08 
0.06 




gure 1: Plot of ||fM<7 = qo,Qo,T) + %±(q = q ,qo,T)\\ where q = (0.01,0) 
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(a) Plot of the solution to the first equa- 
tion defined by Eq. 13.331 



(b) Plot of the solution to the second 
equation defined by Eq. 13.331 




(c) Superposition of the two sets of solu- 
tions 



Figure 2: Periodic orbits for the nonlinear motion about a Libration point 



(a) Plot of the solution to the first equa- 
tion denned by Eq. l-HSl for t = 3.033 + 
0.0005n n e {1 • • • 10} 



(b) Plot of the solution to the second 
equation defined by Eq. 13.331 for t = 
3.033 + 0.0005n n e {1 • • • 10} 




(c) Superposition of the two sets of solutions for t = 3.033 + 0.0005ra n £ 
{1---10} 



Figure 3: Periodic orbits for the nonlinear motion about a Libration point 
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Figure 4: ^/|Ap | 2 + |Ap| 2 q for t e [6days, 35days] 
lunit < > 432m. s^ 1 
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Value of N for t= 59days 

64 days 
70 days 
76 days 
82 days 
88 days 




Figure 6: ^Apd 2 + \&p\ 2 q for t e [59days, 88days] 
lunit < — > 432m. s^ 1 
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